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We study the electronic structures of ABA (Bernal) stacked multilayer graphenes in uniform 
perpendicular electric field, and show that the interplay of the trigonal warping and the potential 
asymmetry gives rise to a number of emergent Dirac cones nearly touching at zero energy. The 
band velocity and the energy region (typically a few tens of meV) of these gate-induced Dirac cones 
are tunable with the external electric field. In ABA trilayer graphene, in particular, applying an 
electric field induces a non-trivial valley Hall state, where the energy gap at the Dirac point is filled 
by chiral edge modes which propagate in opposite directions between two valleys. In four-layer 
graphene, in contrast, the valley Hall conductivity is zero and there are no edge modes filling in the 
gap. A nontrivial valley Hall state generally occurs in asymmetric odd layer graphenes and this is 
closely related to a hidden chiral symmetry which exists only in odd layer graphenes. 



I. INTRODUCTION 

Graphene is characterized with Dirac quasiparticles 
in the low-energy region, ^""^ which give rise to anoma- 
lous physical properties due to the linear dispersion and 
non-trivial Berry phase. ^^^^ There are growing interests 
in multilayer variants of graphene such as bilayer^^-ie 
and trilayer, ^^"21 which also support chiral quasiparti- 
cles. Bilayer graphene has parabolic valence and con- 
duction bands touching each other at Dirac point22"24^ 
while ABA(Bernal)-stacked trilayer graphene comes with 
a superposition of effective monolayer-like and bilayer- 
like bands. 2'^'25-33_ ^-^p of these, the trigonal-warping 
deformation of the energy band, which is intrinsic to 
graphite-based systems^''"'^^ , gives rise to small Dirac 
cones near Dirac point in these multilayers. The Lif- 
shitz transition, in which the Fermi circle breaks up into 

separate parts, takes place at a small energy scale around 
a few meV.22.24,29-31 

On the other hand, it is possible to modify the band 
structure of multilayer graphenes by applying an elec- 
tric field perpendicular to the layer, using external gate 
electrodes attached to the graphene sample. For bi- 
layer graphene, a perpendicular electric field opens a 
band gap, 16,22, 23, 27,37-40 ^j^jjg ABA-stacked iV-layer 

graphenes with > 3, in contrast, some electronic states 
remain around zero energy even in a perpendicular elec- 
tric field, causing an increase of the conductivity near the 
charge neutral point. 

In this paper, we closely study the band structures 
of ABA-stacked multilayer graphenes in the presence of 
uniform perpendicular electric field, and find that the in- 
terplay of the trigonal warping and the potential asym- 
metry generally gives rise to a number of additional Dirac 
cones nearly touching at zero energy, as depicted in Fig. 
1(b) for trilayer graphene. For these gate- induced Dirac 
cones, the band velocity and the energy region (i.e., 
the distance between Dirac point to the Lifshitz tran- 
sition point) are tunable with gate bias voltage. The 
energy region is typically a few tens of meV, which is 
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FIG. 1: Band structures of ABA trilayer graphene in a 
perpendicular electric field depicted in 3D plot (left panel) 
and contour plot (right), for the band models (a) only with 
u,7i,U3 terms and (b) with full parameters. The interlayer 
potential asymmetry is set to be A = 200 meV. 



by order of magnitude greater than in the original non- 
biased multilayer graphene. In a magnetic field, there 
arise triply-degenerate Landau levels originating from off- 
center gate-induced Dirac cones, with wide energy spac- 
ings due to the linear dispersion. 

The gate-induced Dirac cones are generally gapped 
at Dirac point by symmetry-breaking terms. When the 
Fermi energy is in the gap, the system is in a topologi- 
cally non-trivial valley Hall state, where electrons at K+ 
and valleys carry opposite Hall conductivities'^2,43^ 
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A manifestation of the valley Hall state is the emergence 
of chiral edge modes at a zigzag interface which trans- 
ports valley pseudo-spins in an analogous way to the spin 
Hall effect"'^'''^ . The valley Hall state and the helical edge 
modes were previously studied for gapped monolayer and 
bilayer graphenes,^^""^^ and also for ABC (rhombohedral) 
stacked trilayer graphenc.^^ We study the edge states 
a semi-infinite zigzag ribbon of asymmetric ABA multi- 
layer graphcnes, and relate the number of the edge modes 
to the valley Hall conductivity which is a bulk property. 
In trilayer graphene, in particular, we find that non-zero 
valley Hall state is realized in a small external electric 
field, and moreover, a topological transition takes place 
at a certain higher electric field, which is accompanied by 
a change of the number of edge channels inside the bulk 
gap. In four-layer graphene, in contrast, the valley Hall 
conductivity is always zero and there are no edge modes 
filling the energy gap. We show that the nontrivial val- 
ley Hall state generally occurs in asymmetric odd layer 
graphenes, and this is deeply indebted to an approximate 
chiral symmetry peculiar to odd layer graphenes. 

Paper is organized as follows. We briefly introduce the 
effective mass model for graphite in Sec. II, and argue 
the trilayer graphene in Sec. HI in terms of the gate- 
induced Dirac cones, the chiral symmetry, the Landau 
level structure and and the edge states. In Sec. IV, we 
study the four-layer graphene as an example of even-layer 
cases without the chiral symmetry. In Sec. V, we argue 
the chiral symmetry in general odd-layer graphenes, gen- 
eralizing the trilayer's argument. The conclusion is given 
in Sec. VI. 



II. MODEL 



We describe the electronic properties of ABA-stacked 
multilayer graphene using Slonczewski-Weiss-McClure 
model^^"'^^ with hopping parameters described in Fig. 2. 
The model includes intralayer coupling 70, nearest in- 
terlayer couplings 71, 73 and 74, next-nearest layer cou- 
plings 72 and 75, and on-site energy asymmetry A', which 
are estimated in bulk graphite as'^^: 70 = 3.16cV, 71 = 
0.39eV, 73 = 0.32eV, 72 = -0.020eV, 75 = 0.038eV, 
74 = 0.044eV, A' = 0.050eV. A' is the energy difference 
between the sites which are involved in the coupling 71, 
and the sites which are not. 



We consider ABA-stacked A''-layer graphene, where 
\Aj) and \Bj) represent Bloch functions at the point, 
corresponding to the A and B sublattices of layer j, 
respectively. If the basis is arranged as 
1^2), 1^2); ■ ■ ■ ; \An), I-Bjv), the Hamiltonian in the vicin- 




FIG. 2: Lattice structure of ABA trilayer graphene with tight- 
binding hopping parameters. The bottom figure is a top- view, 
and the right is a schematic diagram of the lattice structure. 
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where tt = ^tt^; -|- itt^, tt = p + eA with A being the 
vector potential arising from the applied magnetic field, 
and ^ = ±1 are the valley indeces for K±. The param- 
eter V — '\/3a7o/(2?i) is the band velocity for monolayer 
graphene, and W3 = •\/3i73/ (2?i) a velocity related to the 
band parameter 73, where a w 0.246 nm is the distance 
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between the nearest A sites on the same layer. Uj de- 
scribes the electrostatic potential on j-th layer induced 
by the external electric field, where we assumed a uniform 
potential gradient in the perpendicular direction, This is 
valid in a few-layer graphene with typically <~ 5. For 
thicker multilayers with N >~ 10, the potential drop oc- 
curs within a few layers near the external gate due to the 
screening by the charge carriers in graphene. ^^'^^ 
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III. TRILAYER GRAPHENE 

A. Chiral symmetry and gate-induced Dirac cones 

The Hamiltonian of ABA-trilayer graphene is given 
by Eq. (1) with iV = 3, where the external potential 
is (Ui, 1/2,113) = (A,0, -A). In the absence of A, the 
Hamiltonian can be block diagonalized into monolayer- 
like band and bilayer-like band.^^ In finite A, these sub- 
blocks are hybridized, and it is then useful to arrange the 
basis as 



{(|Ai)-|^3»/%/2, m) + \B3})/V2, 

- |S3))/V2, 1^2), (|Ai) + |^3»/V2},(2) 

where the monolaycr-like band correspond to 1st and 4th 
bases and bilayer-like band to 2nd, 3rd, 5th, and 6th. We 
categorize the first three bases in Eq. (2) as group o, and 
the last three as group •. The Hamiltonian is written in 
this basis as 
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If we keep only relevant band parameters 70, 71, 73 
and the potential A, and neglect remaining parameters, 
the Hamiltonian Eq. (3) possesses the chiral symmetry 
(sublattice symmetry) , in that the diagonal matrix blocks 
Ho and H, all vanish, leaving the off-diagonal blocks D± 
which connects the bases of o to the bases of •. The 
energy spectrum of this simplified Hamiltonian contains 
a single center Dirac cone and six off-center Dirac cones 
at zero energy, as depicted in Fig. 1(a). The robustness 
of gapless spectrum is protected by the chiral symmetry. 
The extra terms with 72 , 75 , V4 and A' in the diagonal 
blocks break the chiral symmetry and open small energy 
gaps at these Dirac cones as in Fig. 1(b). 

The positions of the Dirac points in the chiral Hamil- 
tonian without Ha and iJ, can be found by solving 



At each Dirac point, the degenerate zero-energy bases 
■00 and -0, are derived from the equation D^\ipo) = and 
D^\ip,) — 0, and the effective Dirac Hamiltonian is given 
by 



H, 
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(^o|i?-|^.) 
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(8) 



keeping the lowest order in the momentum shift from the 
Dirac point. By rotating (x, y) coordinate and the spinor 
space by angle 9 at the same time, this is transformed to 



HcS — VxTTxCFx + VyTTyay, 



(9) 



where and ay are Pauli matrices, and {TTx,TTy) is the 
momentum measured from each Dirac point. We find for 
the Dirac point at p = (no need for rotation), 



^0 = (0,1,0) 



^. = (-V27i,A,0)/v/A2+272, 



and 



Av 
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with the valley index ^ = ±. For the off-center Dirac 
points at p = p± , 

00 = {~Ae''',TV2juvp±e^'')/Ci, 
4,, = {~A,TV2v3P±,vp±e~'%/C2, 



Ci = ^A^ + 272 + v^pI,C2 = yjA^ + (2«| + v^)pI, 

(12) 

and the velocities after rotation of (x, y) by 9 become 



2^v{±A'' +^iV3p±)/{CiC2) 
-67iW3up±/(CiC2), 



(13) 



where Vx and Vy correspond to the radial and azimuthal 
directions respectively, with respect to p = 0. and Vy 
of each Dirac points are plotted in Fig. 3(d). The veloci- 
ties are mainly enhanced by applying A, while Vy for p+ 
decreases only slowly. This indicates that the conductiv- 
ity, which is roughly proportional to the square of the 
band velocity^, is enhanced by A, as is consistent with 
the previous transport measurement^* and the theoreti- 
cal cstimation'^^. 
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We have a topological change at gap closing point, A = 
Ac. The Hall conductivity have opposite signs between 
two valleys K± due to the time-reversal symmetry, so 
that the net Hall conductivity is always zero. Neverthe- 
less, the single- valley Hall conductivity is directly related 
to the number of chiral edge modes appearing in zigzag 
edge, as we will see in Sec. HIC. 



TABLE L Sign of chirality, mass and in-gap Hall conductivity 
for gate-induced Dirac cones at valley. 



The chirality for each Dirac cone can be defined by 
i/j, = sgn^v^Vy) , and this coincides with the Berry phase 
in units of tt around the Dirac point. We find Vc = — ^ for 
p = 0, and Vc = foi' P = P±j so that the summation of 
chirality over seven Dirac points in the valley is — ^. 
Since the chirality is a topologically protected number 
as long as the chiral symmetry is present, non-zero total 
chirality indicates that the conduction band and the va- 
lence band inevitably touch at some points in any value 
of A. 

Now we consider all the band parameters in Eq. (3) to 
argue the energy gaps at the Dirac points. The effective 
Hamiltonian of each gate-induced Dirac cone is modified 
to 



(14) 



where the mass m and the energy shift eo are given by 
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The width of the energy gap is 2m. Note that the ap- 
proximation using the gapped Dirac Hamiltonian above 
is valid when the mass gap is smaller than the energy 
region of the gate-induced Dirac cone below the Lifshitz 
transition point. In Fig. 3(c), we show the evaluated mass 
m for the seven Dirac points. We see that the mass for 
p = P- changes from negative to positive when A exceeds 
the critical value. 



Ap « 270 meV, 



(17) 



at which the energy gap closes. 

When the Fermi energy lies in the gap in massive 
Dirac Hamiltonian Eq. (14), the Hall conductivity takes 



[e^ / (2h)]sg\i{vxVym) 



even m zero 



non-zero value (Jxy 

magnetic field. ^^'^^ Considering chirality and mass for 
each Dirac point in Table I, the total Hall conductivity 
summed over the Dirac points at single valley if^ be- 
comes 



--e (A<A,), 
+U (A>A,). 



(18) 



B. Landau level structure 

The Landau levels in the presence of a uniform mag- 
netic field B = (0, 0, B) can be calculated by the Hamil- 
tonian with TT = {\/2'h/t)a) and {y/2Ti/ll)a for K-\. and 
, respectively. Here and a are raising and lowering 
operators, respectively, which operate on the Landau- 
Icvel wave function 0„ as a0„ ~ -y/n(/)„_i, a^^n — 
\Jn -f 10„_|_i, and € = ^jTij {eB) is the magnetic length. 
The Landau level spectrum at B = 0.5T is plotted as 
a function of A in Fig. 3(b). In increasing A, the gate- 
induced Dirac pockets accommodate more and more Lan- 
dau levels as the energy region expands. The level spac- 
ing inside the Dirac pockets is much wider than outer 
region due to the linear dispersion. 

The low-energy Landau levels inside the gate-induced 
Dirac pockets can be approximately described by the 
massive Dirac Hamiltonian, Eq. (14). There the Landau 
level spectrum is explicitly written as 



eo = -sgn{vxVy)m, 
e±n = ±^/A%n + n 



where 



\J2\vxVy 



heB. 
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We note that the Landau level of n = is sensitive to the 
chirality and the sign of the mass: it appears at the top 
of the valence band and the bottom of the conduction 
band when sgn{vxVym) > and < 0, respectively. The 
two cases correspond to the mid-gap values of the Hall 
conductivity, axy = e^/(2/i) and — e^/(2/i), respectively. 
Also, since the chirality sgii{vxVy) is opposite between 
and K_ valleys, the n = Landau level split in 
valleys, while all others are valley degenerate. Besides, 
we have additional triple-fold degeneracies for each of p^ 
and p- . 

In Fig. 3(b), we actually see that the levels n = level 
is non-degenerate in valleys, appearing at either of the 
band edges of each massive Dirac band. At A = Ac, we 
observe that n = levels of K± valleys cross at the charge 
neutral point, in accordance with the topological change 
of crfjf from {-5/2)^e^/h to {+l/2)^e^/h. The levels 
71 > are almost valley-degenerate while tiny splitting is 
due to the deviation from the massive Dirac Hamiltonian 
Eq. (14). 

When we drop the band parameters other than 70,71 
and 73; the Hamiltonian becomes chiral symmetric and 




FIG. 3: (a) Band structure of gated ABA trilayer with various gate voltage A. (b)Gate bias dependence of Landau level energy 
in _B = 0.5T. Landau levels in K+ and K~ valleys are plotted as solid (blue) and dashed (red) lines, respectively, (c) Effective 
mass m and (d) band velocities Vx, Vy for the effective Dirac cones at p = 0,p+,p_, plotted against A. 



the zcro-th Landau level Eq of each Dirac cone comes 
exactly to zero energy as a chiral zero mode, of which 
wavefunction has amplitude only on o and • for Vc = T^i 
respectively. The index theorem then states that the dif- 
ference between the number of the zero-modes belonging 
to o (n+) and those to • (f^-), is defined as chiral in- 
dex, which coincides with the gauge flux penetrating the 
system. The chiral index in the present case is shown 
to be n+ — n_ = where $ = eBS/h is the magnetic 
flux penetrating the system area S. This is, in units of 
$, coherent with a summation of —Vc in each single val- 
ley. As in the conventional Dirac Hamiltonian,^"^"^^ the 
chiral index can be related to the geometric curvature of 
the gauge field, and the above relation between the chiral 
index and total magnetic flux stands in non-uniform mag- 
netic field as well. The detailed argument is presented in 
Appendix. A. 



C. Edge modes 

Non-trivial Hall conductivity in single valley indicates 
an existence of chiral edge modes localized at the inter- 
face, as long as the valley mixing is not present. There the 
number of emergent edge modes are directly related to 
the Hall conductivity, so that chiral edge modes as many 



as the number of cTxy^ should countcrflow in opposite di- 
rections between K-^- and , as is analogous to the spin 
Hall insulator.^"* Here we numerically examined the edge 
modes in the asymmetric ABA trilayer graphene with 
zigzag interfaces, in which the valley mixing is absent. 
We consider a semi-infinite system with a zigzag bound- 
ary along X direction as shown in Fig. 4(a), where Px is a 
good quantum number. The energy of edge modes in the 
bulk gap can be obtained by searching for the evanescent 
modes satisfying a boundary condition at the interface. 
The method is detailed in Appendix B. 

First we consider the chiral symmetric case neglecting 
72,75,i'4, A'. Fig. 4 illustrates the energy spectrum near 
point at A = 200 meV, where we see that zero-energy 
edge modes appear between some of gate-induced Dirac 
points. The number of zero energy edge modes are closely 
related to the chirality of each Dirac cone.^^"^^ When we 
regard a two-dimensional periodic system on xy-plane as 
a one-dimensional system with fixed px as a parameter, 
we can define a winding number "f{px) by integrating the 
Berry phase change all the way along py on the Brillouin 
zone at the fixed Px ■ When we set the boundary perpen- 
dicular to y axis (i.e., Px is still a good quantum number), 
the number of zero-energy edge modes appearing at the 
boundary coincides with j(j>x) except a constant. In 
the present system, this bulk-edge relationship can be 
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FIG. 4: (a) Atomic structure of a zigzag ribbon of ABA tri- 
layer graphene. (b) Band structure around K+ valley of a 
zigzag ribbon of ABA trilayer graphene only including 70,71 
and 73, with the asymmetric potential A = 200 meV. Bulk, 
left-edge and right-edge states are plotted with solid lines, 
red dots and blue dots, respectively. Numbers indicate the 
degeneracies of the edge state bands for a single side (the 
same for the left and the right edges), (c) Contour plot of 
the bulk band structure with white dots and labels denoting 
the position and the chirality of gate- induced Dirac points, 
respectively. One-dimensional winding number 'y{kx) is indi- 
cated between dotted lines penetrating the Dirac points. 





FIG. 5: Band structure around valley of zigzag ribbon of 
ABA trilayer graphene with the parameters fully included, in 
the asymmetric potential (a) A = 200 meV and (b) A = 400 
meV. Bulk, left-edge and right-edge states are plotted with 
solid lines, red dots and blue dots, respectively. 



clearly seen in Fig. 4. 

Other hopping terms breaking the chiral symmetry 
give rise to mass to the Dirac points, and relative signs 
between these masses determines the connection of the 
chiral edge modes between different Dirac points. In 
Fig. 5(b), we plot the band structure near A'+ including 
full band parameters at A = 200 meV. We see that the 
left and right edge modes stick to cither of the top or 
bottom of gapped Dirac cone depending on the sign of 
the mass. At the charge neutral point, we have three set 
of chiral edge channels crossing the Fermi energy, which 
all circulate in clockwise direction when viewed from +z 
direction. We also observe two edge states extending out 
of the plot and leading to the other valley A'_. In 
we have the exactly same spectrum with px inverted to 

The correspondence to the single-valley Hall conduc- 
tivity can be understood in a similar argument to that 
for integer quantum Hall effect. Let us consider a cylin- 
drical system which is closed in x with circumference 
while finite in the axial direction y with —Ly/2 < y < 
Ly/2 bound by the zigzag edges. When we adiabatically 
turn on a magnetic flux quantum h/e penetrating into 
the cyHnder (inducing an electric field along —x direc- 
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tion), every state at is shifted to Px + 2TTh/Lx. The 
single- valley Hall conductivity of then coincides with 
the total move of electrons in — y direction through 
this adiabatic process. At the Fermi energy, an electron 
moves from the left edge to the right for each of three 
pairs of counter-propagating channels, contributing to 
c'xy' = —Se^/h. A charge transfer also occurs below the 
Fermi energy, where a bulk state ((y) = 0) is pumped 
to an edge state ((y) = —Ly/2) in the left-edge chan- 
nel going out of the valley. This yields a contribution to 
crfj+ = {+l/2)e'^/h, which adds up to erf/ {-5/2)e'^/h 
all together. 

Fig. 5(c) plots the band structure at larger bias, A = 
400 meV after the topological transition at Ac. We find 
that the connection of the edge modes changes at p = p_ , 
leaving only one clockwise and one anti-clockwise chiral 
edge modes crossing at the Fermi energy. The Hall con- 
ductivity from those two exactly cancel out, while we 
have the same contribution from the edge channel be- 
low the Fermi energy, giving the total Hall conductivity 
o'xy^ = {+l/2)e'^/h. This again coincides with bulk val- 
ley Hall conductivity estimated from the mass and chiral- 
ity. Conversely, the single valley Hall conductivity gives 
the number of the counter edge modes crossing at the 
Fermi energy, when we appropriately exclude the half in- 
teger contribution from the edge modes connecting K±. 



IV. FOUR-LAYER GRAPHENE 

Unlike trilayer, a four-layer graphene with interlayer 
asymmetry does not possess the chiral symmetry even 
in the approximate model, and the band gap always 
open at the zero energy. The effective low-energy 
Hamiltonian can be obtained by excluding the high- 
energy bonding states at \e\ ~ 0(71) as H^s = Hn — 
H12H22 H2i,^'^ where Hn and H22 represent diagonal 
blocks of the original Hamiltonian for low-energy bases 
spanned by (Al, S2, A3, i34), and for high-energy bases 
by A2, 53, A4), respectively, and H12 and H21 are 
off-diagonal blocks connecting them. This is explictly 
written in basis (Al, B2, A3, 54) as 
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TT TT 
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where we neglected the band parameters other than 70 , 71 
and 73. The approximation is valid only when vp^ A <C 

71- 

If we even neglect term, the low-energy energy band 
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FIG. 6: Band structures of ABA-stacked fourlayer graphene 
with the interlayer asymmetric potential A = 100 meV, de- 
picted in 3D plot (left panel) and contour plot (right). 



is rotationally symmetric around K± points and its dis- 
persion relation is given by 



e{p) = ±- V + 5A2 - 2^/b^^+2^^^^+I^^ (22) 



with e = w p /71. The band gap appears between e = 
±(A/2)(— 7+ 16/a/5)^^^, corresponding to an off-center 
momentum po = {—2 + 6/a/5)"'"^^ VTiS/w. 

When we resume term and other parameters, six 
off-center pockets emerge at momentum p_|_ and p_ near 
Po, each of which are arranged in 120 degrees symmetry 
as illustrated in Fig. 6, and Fig. 7(a). The pre-existent gap 
never closes during this process. The pockets at p = p+ 
are much deeper than those at p = p_, and the energy 
depth is about 15 meV at A = 100 meV. Fig. 7(b) de- 
scribes an evolution of Landau level energies with increas- 
ing A, where we observe the triply-degenerate Landau 
levels of p+ pockets with wide energy spacing, similarly 
to trilayer graphene. 

We can expand the Hamiltonian with respect to the 
center of each Dirac pocket, and obtain the effective 2x2 
Hamiltonian in the massive Dirac form of Eq. (14). The 
masses and the velocities are shown in Fig. 7(c) and (d), 
respectively, where we set the basis so that the chirality 
becomes -1-1 (i.e., VxVy positive). The masses at p+ and 
p_ are opposite in sign, so that the contributions to the 
Hall conductivity of Dirac cones cancel out in summation 
over a single valley. Therefore, the single valley system is 
a trivial insulator with zero Hall conductivity in contrast 
to trilayer graphene. Accordingly, we observe no Landau 
level crossing at charge neutral point in Fig. 7(b), and 
if we look at the edge states in Fig. 7(e), there are no 
edge channel crossing in the bulk gap, nor the counter- 
propagating flows of valley pseudo-spin. 



V. GENERAL ODD-LAYER GRAPHENES 

The approximate chiral symmetry in the presence of 
the external electric field argued in trilayer actually holds 
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FIG. 7: Plots similar to Fig. 3 for ABA fourlayer graphene. (a) Band structure with various gate voltages A. (b) Gate bias 
dependence of Landau level energy in _B = 0.5T. (c) Effective mass m and (d) band velocities Vx,Vy for the effective Dirac 
cones at p+,p_, plotted against A. (e) Band structure similar to Fig. 5 for zigzag ribbon of ABA fourlayer with A = 200meV. 



in any odd-layer Bernal multilayer graphenes. As we see 
in the following, if we only consider 70,71,73, A terms , 
the Hamiltonian of odd layered graphene is chiral sym- 
metric with a nonzero chiral index in the presence of 
magnetic field, which means that the band is gapless in 
limit of zero magnetic field. 

The Hamiltonian for Bernal stacked A'^-layer graphene 
is decomposed into block diagonal form with effec- 
tive monolayer and bilayer blocks with a unitary 
transformation^^''^^''^^''^'^. Let us define 



fni{j) — Cn 



5m (j) = C„ 



iV+ 1 



[1 - (-!)■'] siuK^j, 



iV+ 1 



[1 + (-1)-'] sin Km j, 



(23) 



where 



2 2(7V + 1)' 

'1/2 (m = 0), 
l/^/2 (m^O), 



(24) 



1,2,...,A^ is the layer index, m is the block index 



given as 

/l,3,5,. 
\0,2,4,. 

Then we take a basis 



,N- 



N 
N 



even, 
odd. 



,odd)\ 



X^cvcn) \ 



N 
N 



(25) 



where X = A or B. A superscript (X, odd/even) indi- 
cates that the wavefunction is nonzero only on the sub- 
lattice X on the layer j = odd/even. 

With the basis set above, the Hamiltonian with A = 
is block diagonalized with m. A block labeled by m 



spanned by {\4>i^'°'^'^^) , 10^'°'^'^^), |(/)m '°™"^), \4>m'' 
is dictated as H{Xm) with A ui — Z COS 

and 



>} 



HiX) 







At;37r\ 







im A71 

A71 VTT 

yAwsTrt VTT y 



rt 



(26) 
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where we only left the terms with 70, 7i, 73, For the case 
of m = 0, only two bases \(j>m''°'^'^^) , \(l)m ) survive due 
to goij) = 0, and the corresponding block Hamiltonian 
is the first 2 by 2 component of Eq.(26). 

Now let us consider an odd-layer graphene with the in- 
terlayer potential asymmetry U{A), i.e., the second term 
of Hamiltonian, Eq. (1). When TV = odd, we can easily 
show that the basis function belonging to the block m is 
either symmetric or antisymmetric with respect to reflec- 
tion in the central layer j = {N + l)/2, and the parity is 
given by (-i)(N-m-i)/2 33 jf ^^^^ jj^^) jj-^ ^^lis 

basis, the matrix clement {(l)^^, ■* (j>,p' = even 

or odd) becomes non-zero only when m' = m + Al + 2{l : 
integer), X = X' , and p = p' , because U is an odd func- 
tion in the reflection, and also diagonal in the original 
site representation. 

Therefore, if we separate the basis functions into two 
groups o, • as: 



(14^'°''^ 1^^^'°""^), l^L""'"™"^)) 

J (o, •, o, •) (to ~ Al) 
^|(.,o,.,o) (to = 4/ + 2). 



(27) 



then the Hamiltonian including only 70, 71 , 73, A has ma- 
trix elements only between o and •, and thus is chiral 
symmetric. Note that the chiral symmetry is not re- 
spected for even-layer graphenes since the basis function 
labeled by m cannot be categorized to either even or odd 
parity, and the interlayer potential term gives rise to ma- 
trix elements connecting blocks to and m -\- Al. 

In the presence of magnetic field, the chiral index, i.e., 
the difference between the number of the zero-modes 
belonging to o (n_)_) and those to • (jt.-), can be eas- 
ily obtained by considering the Landau level spectrum 
with 71 , U3 and A all switched off, since the chiral index 
never changes in such a continuous transformation. The 
bilayer-likc Hamiltonian block, Eq. (26), is then consists 
of two monolayer-type 2x2 diagonal blocks, giving two 
zero-energy Landau levels localized at the second and 
fourth (first and third) elements in the valley ^ = -!-(—). 
Considering the base grouping in Eq. (27), the differ- 
ence in zero energy states in block m is found to be 
ri+ — = ^2^3" for m = 4/ and to — Al + 2, respectively, 
where $ ~ eBS/h is the magnetic flux penetrating the 
system. The monolayer-type block to = lacks the third 
and fourth elements, giving n+ — n_ = — As a result, 
the total chiral index in odd-layer graphene for valley 
is finally given by 



n+-n_ = (-1 + 2 -2 + 2- ••)'?* 
f-^$; A^ = 4/ + l 
|+^$; A^ = 4/ + 3. 



(28) 



This states that at least one Landau level remains at 
zero energy, and thus in zero magnetic field, the con- 
duction and valence band touch at one fc-point at least. 
The sum of Berry phases in a single valley coincides with 



— TT times the chiral index in units of that is ±^7r for 
N = Al-\-l and 4/ -I- 3, respectively. When the Dirac cones 
are gapped by including the additional band parameters, 
the Hall conductivity per single valley must be non-zero, 
because the number of Dirac cones per valley must be 
odd to achieve total Berry phase ±7r. 

The result might seem to contradict with the well- 
known fact that the Berry phase is Nt: in A'^-layer 
graphene at A 0, but this is for a different basis group- 
ing in which A and B sublattices belong to o and •, 
respectively. The Hamiltonian with A = is block diag- 
onalized into monolayer-like and bilayer-like blocks, and 
the total Berry phase varies depending on which bases 
are assigned to o or • in each block. In the presence of 
nonzero A, the grouping of Eq. (27) is the only possible 
way to make the Hamiltonian chiral symmetric. 



VI. CONCLUSION 

In Bernal multilayer graphene more than three layers, 
an interplay of the gate electric field and the trigonal 
warping effect gives rise to emergent Dirac cones in the 
low energy bands, whose band velocity and Lifshitz tran- 
sition energy are tunable by the gate voltage. In trilayer 
graphene, in particular, the low-energy effective theory 
shows that the valley Hall state is realized at the charge 
neutral point, where single valley Hall conductivity is 
quantized at a non-zero half integer. We have investi- 
gated the edge states at the zigzag interface, and demon- 
strated that the number of edge modes is closely related 
to the bulk single valley Hall conductivity. In four-layer 
graphene, gate-induced Dirac cones also appear, though 
the system is a trivial insulator with zero valley Hall con- 
ductivity. The non-trivial valley Hall state is generally 
found in odd-layer graphenes, where the approximate chi- 
ral symmetry is responsible for the emergence of non-zero 
valley Hall conductivity. The gate-induced Dirac cones 
should be experimentally accessible directly by observing 
Landau levels with wide energy spacing^^"^^ and also by 
measuring transport properties through the edge modes 
at a zigzag interface. 
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Appendix A: Index theorem for odd-layer graphenes 

Here we show that the chiral index of general odd- layer 
graphenes can be written in terms of the total gauge flux 
penetrating the system, in a similar way to the argument 
for usual Dirac Hamiltonian. ^^^^^ The chiral symmetric 
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Hamiltonian of odd-layer graphene with intcrlayer asym- 
metric potential is written as 

where D+ = {D-Y is a TV x TV matrix. The chiral sym- 
metry is then expressed by TH + HV = with the chiral 
operator F, 

Y=I^N \ 

\0 -In)' 

where Ijv is TV x TV unit matrix. We can take zero modes 
of H as an eigenvector 'ip± of T with eigenvalue ±1, re- 
spectively. If we write the number of chiral zero modes 
'ip± as n±, the chiral index v is defined as 

v = n+-n^= TrTf{H^/M^), 

where / is a regularization function, which is smooth 
and monotonically decreasing function with /(O) = 1 and 
/(oo) = 0, and M is an ultraviolet cutoff. The action of 
/ to the matrix is defined through its action onto the 
eigenvalues, as seen if we take f{z) = e~^. 

If we take plane waves exp(ifc • a;) as a basis for the 
spatial direction, we have 

= Tiif{D^D+/M^) - fiD+D^/M')] 

= /^'^/ (02tre-'=-[/(i?_i^+/Tl/^) 

-f(D+D^/M^)]e''''' 

where Tr means a trace over all the states while tr is a 
trace over the layer and site indeces or, equivalently, m, 
X and even/odd indices in Eq.(25). 

The operator tt acts on a plane wave to give 
g-ikx^gikx _ Ji]^ ^ Since D± contains at most first 
order derivative terms, we can write 

^-^kxJ^_^^kx ^ Dl\k\+D^, 

where Do is a c-number matrix and of zero-th order in 
|fc|, but may depend on a polar angle of k. Since Dq 
is a square matrix, it can be written in a singular value 
decomposition 

Dq = UY.V'^ 

with unitary matrices U, V and a diagonal matrix S = 
diag(yA;) (A, > 0). 

The action of D^D+ on a plane wave is then described 

as 

e-'^'-" D^D±e^'' = G±k^ + F±, 



where 

G+ = DlDa = FS^Ft 
G- =DqDI = UY?U\ 

are c-number matrices of zero-th order in |fc|. F± are 
matrices including the operator tt, and up to first order 
of |fc|. 

Having in mind that the ultraviolet behavior (Tivk ^ 
M) is important for the contribution of D^D^ term to 
the chiral index, we obtain 

tr{e-''=^/(i:>_D+/Tl/2)e'''^} 

= tr{/(fc2EVM^) + f'{k^Y?/M^){V^F+V)/M^} 

+0{M-*) 

= Y.{f{k'K/M') + r{k^K/M^){V^F+V)u/M^)} 

i 

With a similar form for Dj^-D^ term, the chiral index 
with Tl/ — > oo reduces to 

x(ytF+F- C/^F_C/)ri/Af2. (A2) 

As argued in Sec. V, the Hamiltonian of multilayer 
graphene with A = is decomposed into block diag- 
onal form with effective monolayer and bilayer blocks 
labeled by m, and when the number of layers is odd, 
A always enters in the off-diagonal blocks. Then G± is 
block-diagonal because it is independent of A, and so U 
and V are. Therefore, although block off-diagonal terms 
in F± arise from A, they do not contribute to the sum of 
(y^F+V — U'^F^U)ii, and thus we only have to compute 
the chiral index for each block separately and add them 
up to obtain the overall chiral index. 

Monolayer block (m = 0). Using the commutation 
relation [7r,7r'''] = —^2heB, 

D^D+/v^ = tt^'tt ^Tv'^ +^neB 

D+D_/v^ = TTTrt = TT^ - £,heB 

G± = h^v^ 

F± = v^ijr- + 2hk ■ tt) ± ^v^heB. 

From Eq.(A2), 

= y d'x J -f^f'{h\^k'/M^)i2^v'heB)/M^ 

where = / cPxeB/h is the magnetic flux penetrating 
the system. When the magnetic field is uniform, this is 
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the Landau level degeneracy of n = Landau level, and 
its sign reflects that the level is assigned to ■0_ (■)/'+) at 
K+{K-) valley. 

Bilayer block (m ^ 0). For simplicity, we set v ~ 1 
and A = 1, and compensate it by redefining M/v — )■ 
MjXv^/v 113 and X-ji/v — > 71. If we rewrite 
the block Hamiltonian Eq. (26) in an order of bases 

(^,odd), I , (A, even), , ,(B,oven), , ,(B,odd), 



H = 




f V3TT TT^\ 

7rt 71 

W37r^ TT 

V ^ 71 0/ 



we have 



■n-'^{l + vl) 7i7r'f+V37r 



t^2 



7i7r + Vsin^) 



7r-(l + u|) 7i7r + i;3(7rT) 



^7i7r' 



V3TT 




The action on the plane wave basis is then given by 



ikx jj jj ^ikx _ ri ;.2 



F 



G+k^ 



F + F 
hPf -F 



1 + wae^' 



^3 



1 



^heB 



for each block (Eq.(27)), we find that the chiral index for 
biased N (odd) layered graphene is given by 

n+-n_ = (-1 + 2-2 + 2....)^$ 

This exactly coincides with Eq. (28), while the present 
argument is more general and valid for non-uniform mag- 
netic field B{x,y). 



Appendix B: Derivation of the edge modes from the 
effective mass Hamiltonian 



When the Hamiltonian is linear in fc, it is possible to 
obtain the edge state energies in the bulk gap, by search- 
ing for the evanescent modes satisfying a boundary con- 
dition at the interface. We consider a 2N x 2N Hamil- 
tonian matrix H{kx, ky) with k = — iV, and assume it is 
linear in fc. It is expressed as 



where is a matrix including tt, of which expression (not 
presented) is not important in the following argument. 

With a relation tr/(A) = tr/(A"^), an explicit calcula- 
tion of Eq.(A2) shows that 

X [1/t {^F + F)V- {V^)\f'^ - F)yT]^^/Af 2 
= jd'^j ^^Y.nk'KlM^){2V^FV)u/M- 

^ ^ i 
i 

= J d^xtr{G^^F) ^ -2£_'i>. 

Thus, the chiral index of bilayer block is twice of that 
of monolayer block, and inclusion of V3 does not affect 
the result. This is naturally expected because the chiral 
index is a topological number and never changes in a 
continuous deformation. 

If we combine these results for monolayer block and 
bilayer blocks noting that the orders of the chiral bases 



H = Aky + B{k,), 



(Bl) 



where A and B is 2N x 2N matrices and A is inde- 
pendent of kx- We assume the system is periodic in x 
and replace k^ with its eigenvalue k^. H is regarded as 
one-dimcnsional Hamiltonian with a parameter k^. The 
Schrodingcr equation, HF{y) ~ £F{y) is transformed to 



iA- 



-B{kx) + £\F{y) = M{kx,e)F{y), 

(B2) 

with M = iA-^[-B{kx) + e]. Let and /^"H" = 
1 , • • • 2 A^) the eigenvalues and right eigenvectors of the 
matrix M(kx,£), The corresponding wave function be- 
comes 



F(")(2/) = cxp(KMy)/ 



(a) 



(B3) 



Generally k*-"^ is a complex number, and the state 
is a bulk mode when k*-"' is pure imaginary, and an 
evanescent mode otherwise. When the bulk spectrum of 
H{kx, ky) is fully gapped at particular k^, and e is inside 
the gap, we have 27V evanescent modes with N modes 
of Rck'") > and another N modes of ReK*^"^ < 0. 
When we consider the half-infinite system in the region 
y > 0, an edge state localized near y = 0, if it exists, 
should be written as a linear combination of the states 
with RcK^"' < 0. When the indeces a = 1, 2, • • • , TV are 
assigned to the modes of Re k*^") < 0, it is written as 



N 



F{y) = E ^^"^ exp(«:(")y)/(" 



(B4) 



a = l 



The boundary condition at y = for the wavcfunction F 
is composed of linear equations with respect to 2iV- 
dimensional vector F{Q). This is written as DF{Q) = 0, 
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with D being a N x 2N constant matrix. For the trilayer 
graphene with zigzag edge, for example, the boundary 
condition at y = is i^i(O) = ^3(0) = ^5(0) ^ 0, so that 
D becomes 



D 



1 0' 
1 
,000010, 



(B5) 



By using Eq. (B4), the boundary condition is written 



as 



X 



C(2) 



/o\ 



w 



(B6) 



where X = X{kx,e) is a TV x TV matrix defined by 



x = D{r'>,f 



f{1) f(2) 



(W) 



(B7) 



The equation has a non-trivial solution when 
dctX{kx,e) = 0. The edge mode energy can be 
found by tracing detX{kx,e) throughout the energy 
gap, for each fixed k^- When more than two edge 
modes are degenerate at the energy e, the number of 
degeneracy is foimd by — Ra.nkX{kx,e). 
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